Time distribution and loss of scaling in granular flow* 
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Two cellular automata models with directed mass flow and internal time scales are studied by 
numerical simulations. Relaxation rules are a combination of probabilistic critical height (probability 
of toppling p) and deterministic critical slope processes with internal correlation time t c equal to the 
avalanche lifetime, in Model A, and t c = 1, in Model B. In both cases nonuniversal scaling properties 
of avalanche distributions are found for p > p* , where p* is related to directed percolation threshold 
in d = 3. Distributions of avalanche durations for p > p* are studied in detail, exhibiting multifractal 
scaling behavior in model A, and finite size scaling behavior in model B, and scaling exponents are 
determined as a function of p. At p = p* a phase transition to noncritical steady state occurs. 
Due to difference in the relaxation mechanisms, avalanche statistics at p* approaches the parity 
conserving universality class in Model A, and the mean-field universality class in Model B. We also 
estimate roughness exponent at the transition. 
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I. INTRODUCTION 



Dynamics of granular materials represents an impor- 
tant practical and theoretical problem. A new theoretical 
approach to the problem of driven granular flow has been 
initiated in the past few years jjj, which is motivated 
by the observed scaling behavior both in the laboratory 
granular piles and in natural landslides |^]-[| . It has been 
recognized that the collective dynamics of grains may 
lead to a self-organized critical (SOC) states Q, charac- 
terized by scaling properties of sandslides (avalanches). 
Moreover, dynamics may depend on various parameters, 
such as dimension and mass of individual grains and qual- 
ity of their contact surfaces, and on the external condi- 
tions. By varying some of these parameters in a con- 
trolled manner, steady states with different characteris- 
tics are reached, and a phase transition to a steady state 
with no long-range correlations occurs when a parameter 
is varied through certain critical value |J. 

Various cellular automata models have been intro- 
duced so far to mimic stochastic variations in the con- 
ditions of toppling [p|-|l4||. One-dimensional rice-pile au- 
tomata with stochastic critical slope rules have been use- 
ful in understanding transport properties of rice piles Q . 
Relaxation rules in these models are a kind of branching 
processes with internal stochasticity. In two dimensions, 
two models studied in Refs. (l^jlj] utilize mixed dynamic 
critical slope (CS) and critical height (CH) rules, moti- 
vated by the observed nonuniversality of the emergent 
critical states in natural landslides (for a recent review 
see Refs. HQ). 

In the present work we extend the study of the mod- 
els of Refs. [GdlOl, which we term model A and B, re- 



spectively. In these models stochastic toppling by the 
CH mechanism is controlled by an external parameter — 
probability of toppling p, which can be attributed to vari- 
ations in the external conditions (e.g., wetting), or by in- 
ternal kinetic friction, determined globally by the quality 
of contact surfaces between grains. In contrast to rice- 
pile models of Refs. P,ff2l, the present models are more 
appropriate for the evolution of landscape, in which a 
variety of erosion mechanisms might be simultaneously 
active. 

Two types of triggering mechanisms of landslides are 
recognized in the literature [Jl6|,[l7]|l4j : (i) soil moisture, 
which is controlled by rainfall and water level, and (ii) 
ground motion, which influences slope variation. The 
local shear stress threshold may depend on both slope 
angle and soil properties. We assume that these trigger- 
ing mechanisms are dynamically correlated. By wetting 
diffusion probability is lowered and grains stick together, 
thus building up local heights. However, when the differ- 
ence between heights at neighboring sites exceeds certain 
limit, the slope mechanism becomes activated. 

A simplified picture of the natural mechanisms of ero- 
sion is taken into account by combined relaxation rules 
for the height transport on a two-dimensional square lat- 
tice, as follows: If at a site local height h(i,j) > h c , 
then the site relaxes with probability p as h(i,j) 
h(i,j) - 2 ; h(i + l,j±) -> h(i + l,j±) + 1 . If for finite p 
some of the local slopes a±(i,j) = h(i,j) — h(i + l,j±) > 
a c , then the site relaxes with probability one by toppling 
one particle along each unstable slope, repeatedly until 
all slopes are reduced below er c . Here (i + l,j±) are po- 
sitions of two downward neighbors of the site (i,j) on a 
square lattice oriented downward. 



*This paper is dedicated to Professor Franz Schwabl on the occasion of his 60th birthday 
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The system is driven by adding grain from the outside 
at a random site along first row and the system is let 
to relax according to the above rules. Another grain is 
added when the relaxation process stops. The internal 
time scale is measured by the number of steps that the 
relaxation process proceeds on the lattice. At time step 
t—1 a site at first row topples after added grain from 
outside. According to the above relaxation rules, one or 
two grains are toppled from that site, which then appear 
at one or two downward neighboring sites. Therefore, 
mass flow is only down. However, an instability may 
propagate to both downward and upward neighbors of a 
toppled site (except for the sites on the first row, which 
have no upward neighbors) , thus triggering four new sites 
as candidates for toppling per each just toppled site. At 
one time step we update in parallel all candidates for 
toppling. This comprises the usual definition of the time 
step in cellular automata models. 

Since the system builds up unstable sites (with respect 
to probabilistic CH rule), the above dynamic rules need 
to be supplemented by an additional rule, which makes 
the difference between two models. In Model A, all sites 
that are visited by an avalanche at least once are consid- 
ered as candidates for toppling until the whole instabil- 
ity dies off. In this way a propagating instability has an 
internal correlation time t c which is determined by the 
dynamics itself. In Model B, we set t c = 1. Therefore, 
only sites which are in the neighborhood of active sites 
at time t may be candidates for toppling in the next time 
step t+1. It should be stressed that, since an avalanche 
is extended object, in both models there are many sites 
which topple simultaneously and which are not neighbors 
in space. In model B next toppled sites are neighbors 
only on time scale but not in space, whereas in model A 
next toppled sites are not necessarily neighbors neither on 
temporal nor spatial scale. However, all toppled sites are 
connected within affected area in space-time dimensions. 
In both models particles are added from the outside only 
on a random site at the first row and leave the system 
at lower (open) boundary. The mass transport is unidi- 
rectional (down). However, since the above rules allow 
an instability to propagate both forward and backward 
on a 2-dimensional lattice, and evolve on an internal time 
scale, both models are essentially 2 + 1-dimensional, with 
extra dimension representing the internal time scale. Dif- 
ferences in the additional relaxation rule lead to different 
emergent critical states, as explained below. 

In Fig. 1 two examples of large avalanches in model A 
(below) and model B (top) are shown for values of the 
control parameter p = p* at the edge of the scaling region 
(p* ~ 0.4, see Sec. IV for discussion). In both models 
multiplicity of topplings at a site (larger number of top- 
plings is marked by darker gray tone) , is induced by the 
instability propagating back and forth due to nonlocal 
relaxation rules. In model A number of candidates for 
toppling at each time step is larger compared to model 



B, due to internal correlation time typically t c 3> 1, lead- 
ing to more efficient relaxation of unstable sites. On the 
other hand, t c = 1 in model B enables building up nu- 
merous unstable sites (with respect to CH rule) for low 
values of p ~ p* . Therefore huge avalanches with perpen- 
dicular extent comparable to the system size (cf. Fig. 1 
(top)) occur frequently, indicating that the anisotropy of 
the relaxation events vanishes at p*. 

In the limit p — 1 both models reduce to the determin- 
istic directed CH model introduced and solved exactly by 
Dhar and Ramaswamy in Ref. [ |T8| . In this limit slopes 
are restricted to cr±(i,j) < 1, and thus CS rule remains 
inactive. 

II. MODEL A: MULTIFRACTAL SCALING 
BEHAVIOR OF LANDSLIDES 

Correlation times t c > 1 in model A are motivated 
by varying toppling conditions after an avalanche com- 
menced, which represents a natural choice in the case 
of long relaxation times, such as geological evolution of 
landslides. It has been shown that this type of temporal 
disorder is a relevant perturbation both for the evolu- 
tion of landslides [B and for directed percolation pro- 
cesses [jl9). In this model each site develops an individ- 
ual time scale of activity, which then contributes to the 
whole event (avalanche) . As a consequence, the distribu- 
tion of avalanche durations Pa(*, L) in the scaling region 
exhibits multifractal scaling properties when the system 
size L is varied, according to the expression: 

P A (t, L) = (L/L o ) 0t(Qt) ; a t = fin ~) / (hi A 

(1) 

In Fig. 2 the distribution of duration of avalanches is 
shown for p = 0.7 and various lattice sizes. In the in- 
set the spectral function 4>t(o-t) vs. ctt is determined by 
the scaling plot according to Eq. (0) , with to = 1/4 and 
Lq = 1/4. The integrated distribution of durations ex- 
hibits a power-law behavior as P(t) ~ t~( Tt ~ 1 ) in the 
entire region p* < p < 1, with the p-dependent expo- 
nent = n — 1, which is shown in the inset to Fig. 
3. Similar nonuniversality with decreasing scaling expo- 
nents with the parameter p are found for the distribu- 
tions of size D(s) ~ s^' 7 " 3-1 ), and mass of avalanches 
D(n) ~ n^ 1 "™ -1 ) (sec Ref. 14 1 for detailed analysis). 
Slopes of various curves in the main Fig. 3 determine the 
dynamic exponent z(p), which is also shown in the inset 
to Fig. 3. For values of the control parameter p below 
a critical value p* rs 0.4 (see below) the critical steady 
states are no longer accessible by the dynamics. 
III. MODEL B: NONUNIVERSAL SCALING IN 
GRANULAR PILES 

For finite correlation times, i.e., by setting t c = 1, 
avalanches have in the average a reduced number of top- 
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plings per site, compared with model A for the same 
value of the control parameter p. This leads to a smaller 
incidence of large avalanches, and thus to increase of the 
scaling exponents with decreasing probability of toppling 
p. In Fig. 4 the probability distribution of avalanche du- 
rations is shown for few values of p in the scaling region. 
On the other hand, for short correlation times the balance 
between the CS and CH toppling mechanisms is altered: 
By lowering p the system builds up heights faster than in 
the case of model A, and thus the CS mechanism becomes 
more effective, and eventually prevails at the boundary of 
the scaling region at p*. We find numerically that scaling 
behavior is lost at p* < 0.5 JT3|. The scaling behavior for 
p* < p < 1 is characterized by nonuniversalp— dependent 
scaling exponents (see inset to Fig. 5 and Ref. 

The scaling properties of the distribution of avalanche 
durations are determined by using the following finite- 
size scaling form 

P B (t,L)~L- ez T(tL~ z ) , (2) 

where 9 = r t — 1 as above, and z is the dynamic exponent, 
which also depends on p. The scaling plots of Pe(f, L) 
for various values of p in the scaling region and for three 
system sizes at each value of p, are shown in the inset 
to Fig. 4. Similar scaling properties are found for the 
distributions of size and length of avalanches (see Ref. 
p3[ for detailed discussion) . In addition to the temporal 
distribution discussed above, here we also concentrate on 
the distribution of mass of avalanches, Pb(ti, L), satisfy- 
ing the scaling form Pb{u,L) ~ n Tn ~ 1 Q(nL~ Dn ), where 
mass n of an avalanche is determined as total number 
of grains that slide during one avalanche. In Fig. 5 the 
distribution of mass of avalanches is shown for few dif- 
ferent values of the parameter p in the scaling region. 
In the inset to Fig. 5 we plot the exponents 8(p) and 
T n {p) — 1, for duration and mass of avalanches, respec- 
tively, and the dynamic exponent z(p), and fractal di- 
mension of mass D n (p) against p. For p > 0.5 the fol- 
lowing scaling relations are satisfied (cf. inset to Fig. 
5): (t„ — l)D n = z8 = a, where a = T( — 1 is the ex- 
ponent of length of avalanches, which is determined in 
Ref. . The dynamic exponent z which appears in the 
scaling form (^) can also be determined from slopes of 
the curves < T >i vs. £, similar as we have determined 
it in model A. Obtained values are in a good agreement, 
within numerical error bars, with those obtained from the 
scaling plots in Fig. 4. Values of the exponents at p = 0.4 
arc taken from the straight sections of the lines represent- 
ing distributions of duration and mass for smaller system 
sizes L < 128. As indicated in the inset to Fig. 5, these 
values do not satisfy scaling relations within error bars, 
indicating that p = 0.4 is already beyond the edge of the 
scaling region in model B (see discussion in Sec. IV). 



IV. UNIVERSAL CRITIC ALITY AT THE EDGE 
OF THE SCALING REGION 



When the control parameter p is varied through a 
critical value p* we find that the scaling behavior of 
the avalanche distributions is lost, indicating that self- 
organized critical states are no longer accessible by the 
dynamics. By numerical simulations of various distri- 
butions and applying the appropriate scaling analysis it 
was shown that critical steady states disappear below 
p = 0.4 in model A and below p = 0.5 in model 
B |lj|. Here we argue that dynamic rules with different 
correlation times t c in these models lead to separate pre- 
vailing relaxation mechanisms at the edge of the scaling 
region, which lead to different values of p* and to two 
different universality classes of scaling behavior. In par- 
ticular, in model A we find that the scaling exponents of 
large avalanches 0(p*), r s (p*), etc., are close to the uni- 
versality class of parity-conserving (PC) branching and 
annihilating random walks |2(||2l]], whereas in model B 
the exponents at p* reach the values of the mean-field 
SOC universality class. 

Although the relaxation rules in both model A and 
model B are complex interplay of the probabilistic criti- 
cal height and deterministic critical slope rules, we may 
distinguish two basic type of local branching and annihi- 
lating processes that take part to propagate an avalanche 
in these models. Propagation of an avalanche may stop 
at a site to which one or two particles drop in time step 
£, in the following two cases: (1) One particle drop will 
not continue if the site had height zero, that is "annihila- 
tion" A — > occurs with probability 1 — p, where p is the 
dynamically changing probability that a site has height 
h > 1; (2) When two particles drop to a site at time t, 
the avalanche may not proceed if the diffusion probabil- 
ity p is too low, i.e., A + A — > occurs with probability 
1 — p. Note that since number of particles is conserved 
by the processes in the interior of the pile, "annihila- 
tion" means accumulation of particles at a site, which 
thus will take part in future events, in contrast to the 
case of chemical reactions, where particles are extinct. 
When the conditions for toppling are fulfilled, propaga- 
tion of an avalanche represents a branching process which 
consists of two steps. A toppled site at time t transfers 
two particles forward, however, the instability is trans- 
ferred to its four neighbors, but the site itself can not 
topple in the next time step. Toppling of an isolated 
site away from the open boundaries by the critical height 
(CH) rule makes four neighboring sites as candidates for 
toppling in the next time step, and if these four sites 
topple, they make nine new candidates for toppling etc, 
along the chain 1^4^9^16^ 25, • • •. Since each 
toppled site, both initial and triggered sites, topple by 
two particles, in CH mechanism, this toppling chain rep- 
resents a reaction A — > (m + l)A with odd number of 
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offsprings m = 3, 5, 7, 9, • • • per each initial particle . The 
same conclusion is true for topplings by the critical slope 
(CS) rule with two simultaneously unstable slopes. If, 
however, a site topple by critical slope (CS) rule by drop- 
ping one particle along one unstable slope, it will trigger 
three neighboring sites to topple by CS mechanism, and 
another toppling chain occurs as 1 — > 3 — > 7 — > • • •, i.e., 
m = 2, 4, 6, • • ■ offsprings per initial particle. 

Diffusion limited branching annihilating random walks 
(BARW) have been studied by field-theory methods (for 
a recent review see ]22j and references therein). It has 
been recognized that d c — 2 is the upper critical space di- 
mension, and that BARW with even number of offsprings 
in d = 1 belong to PC universality class, whereas the 
directed percolation (DP) universality class was found 
in the case of odd number of offsprings. Two examples 
of dynamical processes in 1 + 1 dimension belonging to 
PC universality class have been studied numerically in 
Refs. 0] and |p4| . It should be stressed that in contrast 
to BARW and DP processes, the present models A and 
B are dynamical, and thus the propagation rules apply 
statistically and depend on the history of the state of 
the system. Recently an analogy between the directed 
percolation and stochastic dynamical model with critical 
height rules has been discussed in Ref. p5[ . 

In model B, short correlation time makes the relax- 
ation at a site less effective with decreasing p, and thus 
efficient building-up of heights occurs, leading eventually 
to p(p*) ~ 1. (A transverse section of the pile in model B 
at p* is given in Fig. 6 (bottom)). Therefore, when an in- 
stability starts, it may trigger a mixture of branching pro- 
cesses described above, making an instability transferred 
back and forth on two-dimensional lattice and evolving 
in time. Fractal dimension associated with the mass of 
avalanches at the edge of scaling region was found to be 
D n {v*) ~ 2 (see inset to Fig. 5). Thus an avalanche 
appears to be compact in 2-dimensional space and, since 
next toppled sites are neighbors in time (t c = 1), it repre- 
sents a connected object in 2 + 1 dimensions. Starting an 
instability in full lattice, i.e., with no threshold condition, 
will trigger an avalanche which propagates as a directed 
percolation cluster in 3-dimensions, until eventually too 
many sites will have heights zero and the avalanche will 
stop. Thus p* should coincide with the site-directed per- 
colation threshold on simple cubic lattice pf DP = 0.435 
p6| . Mass of avalanche is defined as the number of parti- 
cles that slide during an avalanche, and thus it is equiva- 
lent to number of branchings. Therefore, since for p = p* 
the effective dimension D n (p) reaches the upper critical 
dimension of BARW, we may expect mean-field univer- 
sality class for the scaling behavior of avalanches. Our 
numerical results listed in Table 1 confirm this conclu- 
sion. A schematic phase diagram is shown in Fig. 7. 

The situation is different in model A (cf. Fig. 1), where 
decreasing diffusion probability p an avalanche is either 
extinct quickly (short avalanches), or lives much longer 



(large avalanches) with large separation times jL4j. In 
turn, this leads to the efficient topplings at each affected 
site due to many attempts within correlation time t c 3> 1. 
In the resulting steady state most of the sites have heights 
h < h c (cf. Fig. 6 (top)). Therefore, only toppling by CH 
rule takes place and threshold condition is still active, in 
contrast to model B. In model A a toppled site at 
time t may trigger topplings at time t + 1 at three neigh- 
boring sites, since the site toppled at t — 1 time step will 
not fulfill the threshold condition (h > 2) at time t + 1. 
It turns that among three neighbors less than two sites 
topple in the average, therefore leading to a chain of top- 
pled sites with few branches, which is embedded in 2 + 1 
dimensional space-time. However, affected sites which do 
not topple due to low probability p at first attempt may 
topple in later time steps before the avalanche dies off, 
thus starting a new chain. The avalanche is made of set 
of such chains, and has the fractal dimension D n = 1.48. 
We believe that this effectively low dimensional BARW 
process, although it takes part in 2 + 1 dimensional space- 
time, is the reason for PC universality class in model 
A. Another reason for the occurrence of PC universality 
class in reaction-diffusion processes might be the exis- 
tence of more than one symmetric absorbing states, as 
discussed in 0] and ^7j. The process is reminiscent to 
bond-directed percolation in 3-dimensional simple cubic 
lattice, thus we also expect that p* < pf DP = 0.382 
p6[ . In the phase diagram in Fig. 7 phase boundaries for 
model A (dashed lines) separate reactive phase from the 
critical and noncritical absorbing phases. 

In the phase diagram in Fig. 7 phase boundaries for 
model A (dashed lines) separate nonconducting phase 
from the conducting critical and noncritical phases. In 
model B the noncritical conducting phase exists only 
along the line p = 1 below MF point, and a finite slope 
occurs via a phase transition at SB point On the 

other hand, in model A our results suggest that noncrit- 
ical steady states occupy a finite region close to the right 
corner, and that a finite slope occurs asymptotically at 
p = 0, p = 1. Further analysis is necessary in order 
to find precise location of the PC point and the nature 
of phase transition between critical and noncritical con- 
ducting states. Along the phase boundaries between the 
points PC and DR, and between MF and DR, we have the 
nonuniversal criticality of model A and model B, respec- 
tively, discussed in the present work. The point marked 
by DR at p = 0.5, p — 1 corresponds to the universal 
SOC of Dhar-Ramaswamy model. 

Sets of the exponents for p = 0.4 are summarized in 
Table I. Exponents in the model B at this value of p are 
estimated from the straight sections of lines in the sub- 
critical region for smaller lattice size L — 128. Value of 
the exponent r s is taken from Ref. Bq], For comparison, 
shown are also the numerical values of the exponents for 
PC universality class, from Ref. |2l|], and mean- field self- 
organized criticality exponents, from Ref. [E8[ . Note that 
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our exponent 9 corresponds to the survival probability 
distribution exponent S in Ref. pl[ |, and z = 2v^jv\\ and 
that the scaling relation t s — 1 — 0/(9 + 1) holds. The ex- 
ponents T n for mass of avalanche and roughness exponent 
X are unique for granular piles, and can not be defined in 
models of chemical reactions or damage spreading, con- 
sidered in Rcfs. [ p3|]2~4j | . We estimate roughness exponent 
X from the contour of several transverse sections of the 
pile (two examples are given Fig. 6). For instance, by 
using box counting method we find the fractal dimension 
of the contour curve in model B, as df — 1.179 — 1.183, 
and using the expression x = df ~ 1 leads to the value 
listed in Table 1. 

V. CONCLUSIONS 



the order parameter j3 and correlation length vu should 
belong to two distinct universality classes of the dynamic 
phase transitions in these models. Our results suggest 
that although basic relaxation rules in laboratory granu- 
lar piles and natural landslides might be the same, details 
of actual implementation of these rules such as variation 
of control parameter during the course of an avalanche 
might lead to entirely different critical states. 
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We have shown that sandpile automata with mixed 
relaxation rules of stochastic diffusion and deterministic 
branching processes are capable of describing nonunuiv- 
ersality of the self-organized critical states and a loss of 
scaling at a critical value of the control parameter, in 
a qualitative agreement with observed behavior in natu- 
ral and laboratory granular flow. Differences in the re- 
laxation rules due to internal correlation time lead to 
distinct dynamic critical states. In particular, unlim- 
ited (within lifetime of an avalanche) correlation time t c 
in model A leads to a multifractal scaling behavior and 
scaling exponents of large avalanches decrease with de- 
creasing values of the control parameter p. On the other 
hand, finite correlation time t c = 1 in model B leads to 
increase of the scaling exponents with decreasing p, and 
to finite-size scaling properties of avalanches in the entire 
scaling region p* < p < 1. At the edge of the critical re- 
gion at p* , dominating relaxation mechanisms of modulo 
two conserving branching processes and effectively low 
dimensionality of the processes, lead to criticality in the 
parity conserving universality class in model A. In model 
B building up of a global slope appears to be dominant on 
top of the above branching processes, which thus appear 
to have the effective dimension which exceeds the upper 
critical dimension of BARW, and thus mean-field scal- 
ing exponents. It should be stressed that the numerical 
values of the exponents listed in Table I prove the close- 
ness of these universality classes within numerical error 
bars, which we estimate as 0.03. Value of the exponent 
r n = 1.66 in mean-field models is known only numerically 
p8| , whereas in branching processes, which are equivalent 
to sandpiles with a fixed number of grains per toppled 
site, there is the equality r„ = t s = 3/2. Study of the 
details of the dynamic phase transition in these models, 
e.g., in terms of the order parameter and its fluctuations, 
is left out of the present paper (see, Ref. |f| for appear- 
ance of finite slope at SB point in model B). However, 
due to scaling relations among various exponents at the 
transition, the observed different universality classes of 
avalanche statistics at p* indicate that the exponents of 
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TABLE I. Critical exponents at p — 0.4 in models A and 
B, and for parity-conserving (PC), and mean-field universality 
classes (MF-SOC). 
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FIG. 2. Double-logarithmic plot of the integrated distri- 
bution P(t) vs. t for p=0.7 and for various lattice sizes L=12, 
24, 48, 96, and 192 (left to right) in model A, obtained by 
open boundary conditions. Inset: Multifractal spectral func- 
tion 4>t{ott) vs. a t . (Fig. 4 from Ref. [14]). 
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FIG. 3. Average duration time of avalanches of selected 
length < t >t vs. length L for various values of p=l, 0.7, 
0.6, 0.5, and 0.4 (bottom to top) in model A. Maximal length 
equals the system size L max = 192, except for p —OA, where 
Lmax =128. Inset: Scaling exponents 6 = t% — 1, z, and r n — 1, 
defined in the text, vs. p in the scaling region. 



FIG. 1. Two examples of large avalanches running from 
left to right at p = p*: in model A (below) and in model B 
(top). Multiple topplings up to forth order are marked by 
different degree of gray color. 
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FIG. 4. Double-logarithmic plot of the integrated distri- 
bution P(t) vs. t for various values of the parameter p=l, 0.8, 
0.7, 0.6, and 0.5 (top to bottom) and system size L —192 in 
model B. Inset: Finite size scaling plots for the same values 
of p as in the main Figure (left to right). For each plot three 
different lattice sizes L —48, 96, and 192 are used. Plots for 
different values of p are shifted to the right for easier display. 
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FIG. 5. Double-logarithmic plot of the integrated distri- 
bution of mass of avalanche D(n) vs. mass n for L=192 and 
for p=l, 0.8, 0.7, 0.6, and 0.5 (top to bottom) in model B. 
Inset: Scaling exponents a = Tt — 1, 6 = r% — 1, and r n — 1, 
fractal dimensions 2 and D n , and products D n (r n ~ 1) — z8 
plotted against p in the scaling region (see text). 
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FIG. 6. Transverse section of the pile at distance i = 48 
from the top of pile for p = 0.435 in model B (below), and for 
p = 0.382 in model A (top). 
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FIG. 7. Schematic phase diagram in (p,p) plane. Phase 
boundaries for model A (dashed lines) and for model B (solid 
lines). 
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